- Fix libxc_tools.cpp vtxc: use chr->rho[is][ir] (nspin=4 density)
  instead of libxc-format rho (up/dn spin density)
- Revert libxc_setup.cpp to develop (the accel branch changes
  had map key collision, PBEH param count error, and weakened
  param validation)
- Revert input_parameter.h: restore develop defaults for device,
  sc_* params, xc_exch_ext/xc_corr_ext; keep only gga_grad=3
- Revert docs to develop version

- libxc_pot.cpp: comments for SF gradient decomposition, gga_grad branching
- libxc_tools.cpp: comment for nspin=4 SF dh path with 4-component potential
- xc_pot.cpp: comment for non-libxc SF builtin dispatch
- xc_functional_gga_noncol_sf_builtin.cpp: header comment block with
  Scalmani-Frisch reference, chain-rule formulas, and step-by-step labels

…thm (gga_grad=0)

- Add gga_grad parameter to XC_Functional::v_xc and
  XC_Functional_Libxc::v_xc_libxc instead of reading PARAM.inp.gga_grad
  inside the xc module; update all call sites.
- Pass use_sf/has_mag down to convert_vtxc_v, convert_v_nspin4 and
  convert_v_nspin4_sf; drop PARAM reads from libxc_pot.cpp,
  libxc_tools.cpp, xc_pot.cpp and the SF builtin (net global-dependency
  delta of the PR is now zero).
- gga_grad=0 (new default) selects the original algorithm so existing
  results are unchanged; gga_grad=1 keeps the collinear approximation,
  2/3 select the SF-based paths.
- Revert unrelated default arguments added to libxc_abacus.h interfaces
  and restore the init_func doc comment.
- Forward-declare Charge in xc_functional_gga_noncol_sf_builtin.h.

…_vtxc_v

The gga_grad=2/3 branch of convert_vtxc_v wrote the 4-component dh into
v(2,:) and v(3,:) of a 2-row matrix (out-of-bounds). Map dh back to the
spin-up/down representation with the inverse of the SF decomposition:
v_up -= dh0 + sum_mu dh_mu * m_hat_mu, v_dn -= dh0 - sum_mu dh_mu * m_hat_mu.
The vtxc accumulation over the 4 nspin=4 channels is unchanged.

Methods 2 and 3 share the same chain-rule spin-up/down gradients
grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but differ in the
divergence of h = df/d(grad rho) in the potential:
- gga_grad=2 (projected): v_mu -= m_hat_mu * div((h_up - h_dn)/2),
  dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms;
- gga_grad=3 (full SF): v_mu -= div((h_up - h_dn)/2 * m_hat_mu),
  retaining all cross terms.
Implement the projected variant in cal_dh_sf (LIBXC path) and in
v_xc_ncgga_sf_builtin (built-in path, now also dispatched for
gga_grad=2), and document the formula-level differences in the code
comments and in the gga_grad INPUT documentation.